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LATTICE BOLTZMANN METHOD FOR 3-D FLOWS WITH CURVED BOUNDARY 

RENWEI MEl*, WEI SHYY*, DAZHI YU 1 , AND I.I-SHI UUO$ 


Abstract. In this work, wo investigate two issues that are important to computational efficiency and 
reliability in fluid dynamics applications of the lattice Boltzmann equation (LBE): (1) Computational stabil- 
ity and accuracy of different lattice Boltzmann models and (2) the treatment of the boundary conditions on 
curved solid boundaries and their 3-D implementations. Three athermal 3-D LBE models (D3Q15, D3Q19, 
and D3Q27) are studied and compared in terms of efficiency, accuracy, and robustness. The boundary treat- 
ment recently developed by Filippova and Hanel and Mei et al. in 2-D is extended to and implemented for 
3-D. The convergence, stability, and computational efficiency of the 3-D LBE models with the boundary 
treatment for curved boundaries were tested in simulations of four 3-D flows: (1) Fully developed flows in 
a square duct, (2) flow in a 3-D lid-driven cavity, (3) fully developed flows in a circular pipe, and (4) a 
uniform flow over a sphere. We found that while the fifteen-velocity 3-D (D3Q15) model is more prone to 
numerical instability and the D3Q27 is more computationally intensive, the D3Q19 model provides a balance 
between computational reliability and efficiency. Through numerical simulations, we demonstrated that the 
boundary treatment for 3-D arbitrary curved geometry has second-order accuracy and possesses satisfactory 
stability characteristics. 

Key words, lattice Boltzmann method, boundary condition for curved geometries, accuracy, 3-D flows, 
Navier-Stokes equations 

Subject classification. Fluid Dynamics 
1. Introduction. 

1.1. Basic Notion of the Lattice Boltzmann Equation. In one fashion or another, conventional 
methods of computational fluid dynamics (C'FD) compute pertinent flow fields, such as velocity u and 
pressure p, by numerically solving the Navier-Stokes equations in space x and time t [23, 9, 27]. In contrast, 
various kinetic methods use the transport equation, or the Boltzmann equation in particular, for various 
problems in fluid dynamics. The Boltzmann equation deals with the single particle distribution function 
/(x, £, t), where 4 is the particle velocity, in phase space (x, £) and time t. Recently, the method of the 
lattice Boltzmann equation (LBE) [5, 24, 3, 6] has become an alternative to the conventional CFD methods 
employing Navier-Stokes equations. The theoretical premises of the LBE method are that (1) hydrodynamics 
is insensitive to the details of microscopic physics, and (2) hydrodynamics can be preserved so long as the 
conservation laws and associated symmetries are respected in the microscopic or mesoscopic level. Therefore, 
the computational advantages of the LBE method are attained by drastically reducing the particle velocity 
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space £ to only a very few discrete points without seriously degrading hydrodynamics. This is possible 
because the LBE method rigorously preserves the hydrodynamic moments of the distribution function /, 
such as mass density and momentum fluxes, and the necessary symmetries [13, 14, 1]. 

One popular kinetic model is the Boltzmann equation with the single relaxation time approximation [4]: 

(1.1) a/ +€-v/ = -![/- /< o >], 

where £ is the particle velocity, / (0 ^ is the equilibrium distribution function (the Maxwell-Boltzmann distri- 
bution function), and A is the relaxation time. The mass density p and momentum density pu are the first 
(D + 1) hydrodynamic moments of the distribution function / and /(°), where D is the dimension of velocity 
space. 

To solve for / numerically, Eq. (1.1) is first discretized in the velocity space £ using a finite set of 
velocities {£ a } without affecting the conserved hydrodynamic moments [13, 14, 1], 

(1.2) dt f a + U • V/ a = - 1 [f Q - /<°>] . 

In the above equation, / Q (x, t) = f(x, £ a ,i) and fi°\x, t) = / (0 *(x,£ a > t) are the distribution function 
and the equilibrium distribution function of the a-th discrete velocity respectively. The nine- velocity (or 
9-bit) LBE model on the 2-D square lattice, denoted as D2Q9 model, has been widely used for simulating 2-D 
flows. For 3-D flows, there are several cubic lattice models, such as the fifteen- velocity (D3Q15), nineteen- 
velocity (D3Q19), and twenty-seven-velocity (D3Q27) models [12], which have been used in the literature. 
All three models have a rest particle (with zero velocity) in the discretized velocity set {£ Q }. A minor 
variation of those models is to remove the rest particles from the discrete velocity set; the resulting models 
are known as the D3Q14, D3Q18, and D3Q26 models, respectively. The LBE models with a rest particle 
generally have better computational stability. For athermal fluids, the equilibrium distributions for D2Q9, 
D3Q15, D3Q19, and D3Q27 models are all of the form [14] 

(1.3) fi eq) = w Q p 1 + ^e Q -u+ ^j(e a -«) 2 - T( u . w ) , 

where w a is a weighting factor and e a is a discrete velocity, c = Sx/St is the lattice speed, and 6x and St 
are the lattice constant and the time step, respectively. (The values of the weighting factor w Q for D3Q15, 
D3Q19, and D3Q27 models and the diagrams illustrating the lattice structure for D3Q15 and D3Q19 models 
are given in the Appendix.) It can be shown that fa 9 ^ is in fact a Taylor series expansion of the Maxwellian 
fL° ] [13, 14]. This approximation of by the above fa Q ^ makes the method valid only in the incompressible 
limit u/c -4 0. 

With the velocity space discretized, the hydrodynamic moments of / and are evaluated by the 
following quadrature formulas: 

(1.4a) *=5> = £/£”. 

a a 

(l-4b) pu = '£e a f a = 'Ze°fL° ) - 

a a 

The speed of sound of the above 3-D LBE models is c 3 = cj\/Z and the equation of state is that of an ideal 
gas p = pc 2 s . The viscosity of the fluid is v — Ac^. 

Equation (1.2) is often discretized in space x and time t into 

(1.5) f a {xj + e a St, t + St) + f a (xi , t) = -t [ f a (xi , t) - f l a 0) (x t , <)] , 
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where r = X/St. This is the lattice Boltzmann equation with the Bhatnagar-Gross-Krook (BGK) approx- 
imation [4] and is often referred to as the LBGK model [5, 24]. The viscosity in the NS equation derived 
from Eq. (1.5) is 


(1.6) v = (t - 1/2 )c 2 s St. 

This choice for the viscosity makes the LBGK scheme formally a second order method for solving incom- 
pressible flows [14], The positivity of the viscosity requires that r > 1/2. Equation (1.5) can be solved in 
the following two steps: 

(1.7a) collision step: fa{xi, f) = fa(xi, 0 — ~ |/<*(*tt 0 — *)] • 

(1.7b) streaming step: /«(ii + c a <!M + ^)-/o( I ii *)> 


where f Q denotes the post-collision state of the distribution function. It is noted that the collision step is 
completely local, and the streaming step is uniform and requires little computational effort. Equation (1.7) 
is explicit, easy to implement, and straightforward to parallelize. 

1.2. Boundary Condition on a Solid Surface. To date, most Neumann type boundary conditions 
for a solid boundary used in the LBE method are based upon the bounce-back boundary condition: A 
particle colliding with a stationary wall simply reverses its momentum. Much of the previous work on LBE 
boundary conditions is devoted to the analysis and improvement of the bounce-back boundary condition 
[29. 10, 17, 2, 22, 7, 11, 15, 19, 30]. The bounce-back boundary condition can attain second-order accuracy 
if the boundary is fictitiously placed half-way between two nodes. That is, the second-order accuracy of the 
bounce-back boundary condition can only be achieved when the boundaries are located right in t in' middle 
of two neighboring lattices [A -- 0.5; see Eq. (1.8)]. (Readers are referred to our recent work [20] for a 
summary of the previous work.) This prevents the direct application of the bounce-back tjpe boundary 
conditions to simulate a solid body with smooth curvature. To circumvent this difficulty, Mei k Shvy solved 
Eq. (1.2) in curvilinear coordinates using a finite difference method to solve for f a [21]. One can also use 
body-fitted curvilinear coordinates with interpolation throughout the entire mesh, except at the boundaries 
where the bounce-back boundary condition is used [12]. In more recent works [8, 20], Cartesian coordinates 
are adopted with interpolation used only at the boundaries. These techniques rely on the freedom of using 
interpolation techniques. We used the latter technique in the present work. 

As shown in Fig. 1 for a 2-D projection involving a 3-D body, the streaming step requires the knowledge 
of fa{xb, <), in which e 5 = -e Q , at x b on the solid side in order to compute / a (*/, t) for the lattice node 
located on the fluid side at Xf = x b + e&5t. Defining 


( 1 . 8 ) 


A = 


Ik/ 


X, 


I \xf - aSftll 

as the fraction of an intersected link in the fluid region, it is seen that 0 < A < 1 and the horizontal or 
vertical distance between x b and x w is (1 — A )8x on the cubic lattice. 

Based on the work of Filippova and Hanel [8], hereinafter referred to as FH, Mei et al. [20] proposed 
the following treatment for on curved boundaries, 


(1.9) 


fa(x 6, <) = (!- x)/«(*/. *) + xtiHxb, t) + 2 w a p 


0^-2 e ° 


U, 


with 

( 1 . 10 ) 


3 9 . , 2 3 , 

1 + -e a -u bf + ^(e a -u f ) -—(uf-Uf) 


Ji ^ OlP 
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Fig. 1. 2-D projection of the layout of the regularly spaced lattices and curved wall boundary. The thick curve marks the 
boundary location. The solid circles (•) mark the positions where particle- boundary collision occurs . The empty (o) and shaded 
(*) circles are fluid sites and solid sides, respectively. 


and 




(1.11a) 

(A - 1) 1 

U bf = A Uf+^U W 

and \ — — (2A — 1) 

T 

for A > 1/2, 

(1.11b) 

Ubf = Uff = u/(xf + e & St, t) 

J (2A - 1) 

“ d * = (r-2) 

for A < 1/2. 


It is noted that Eq. (1.11b) for and \ differs from that originally proposed by FH. The choice for Ubf 
given by Eq. (1.11b) improves the computational stability for r < 1 and A < 1/2 [20]. Since Eqs. (1.9) - 
(1.11) are in vector form, they can be directly extended to 3-D flows with curved boundaries. 

1.3. Scope of the Present Work. The present study examines two issues in 3-D incompressible fluid 
dynamics simulations with arbitrary boundaries using the LBE method: i) The performance of various 3-D 
athermal LBE models for viscous flows, and ii) the efficacy and reliability of the extension of the curved 
boundary treatment from 2-D to 3-D flows. We focus on the stability and accuracy of the computation and 
the robustness in handling an arbitrary curved geometry. In Section 2, a modification of the choice of u^f 
and the expression for \ when A > 1/2 is proposed in order to further improve the computational stability 
of the boundary treatment. In Section 3, numerical results for four 3-D steady flows are examined and 
various computational issues are addressed. These four cases are: (i) pressure driven fully developed flow in 
a square duct; (ii) 3-D lid-driven cavity flow; (iii) pressure driven fully developed flow in a circular pipe; and 
(iv) uniform flow* over a sphere. In cases (i) and (iii), the LBE-based numerical solutions can be compared 
with known exact solutions so that the accuracy of the LBE solutions can be determined. The difference in 
these two cases is that A is a constant in the square duct while A varies around the solid boundary in the 
circular pipe. In the lid-driven cavity flow 7 , the singularity at the corners between the moving and stationary 
walls allows for a performance assessment of various LBE schemes. The flow past a sphere is an external 
flow around a 3-D blunt body. In all four cases, detailed assessments are made in terms of error norms and 
velocity profiles. It will be demonstrated that accurate and robust solutions are obtained using the newly 
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proposed boundary conditions along with the selected LBE models. 

2. Modification of the Boundary Condition for A > 1/2. Equations (1.9) (1-11) arc first applied 

to a fully developed pressure driven 2-D channel flow by using the 3-D LBE model D3Q19. At the inlet 
(j = i) and exit (i = N x , in which N t is the number of lattices in the ^-direction) the following zero derivative 

condition is imposed after the collision step, 

(2.1a) /<*(* = E E k) = /o(* = 2, E ^), 

(2.1b) /q(® = /■ k) = f a (i — A'a- — 1, j, k). 


At k = 1 and k = N z , the same is imposed, 


(2.2a) fa(i, j,k=l) = f a {i, j, k = 2), 

(2.2b) U(i, j, k = N z ) = f a (i, j, k = N : - 1). 


The constant pressure gradient Np along the ^-direction is treated as a body force and is included in the 
solution procedure after the collision step and the enforcement of the above zero-derivative conditions as: 


3 dp 

(2.3) fa(Xi, t) = f a (x i, t) - w a -j—e a -x, 

where x is the unit vector along the x-axis. On the solid walls (y = 0 and y = H), Eqs. (1.9) - (1 11) are 


used. The exact solution for the velocity is used as the velocity initial condition. The equilibrium distribution 
fa 9) function based on the exact solution for the velocity profile is used as the initial condition for / „ . The 
pressure gradient is set to g = -1.0 x lO" 6 . All computations are carried out using double precision. 



A 


Fig. 2. Stability boundary of FH’s scheme in a square duct flow for A near 1. Shaded areas under stability boundaries 
are unstable regions. 

It was found that the computations are stable for r close to 0.5 (for example, t = 0.505) as long as A 
is not too close to unity (for example, A < 0.87). When A is equal to 1, stable computation can only be 
carried out for r no smaller than 0.6. Fig. 2 shows the stability bound for the channel flow simulation with a 
system size N x x N y x N. = 5 x 35 x 5, near A = 1. Also shown by the dashed line is the stability-instability 
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boundary for the channel flow simulation using D2Q9 model and with a system size N x x N y = 5 x 35, near 
A = 1. It is clear that similar behavior exists in both 2-D and 3-D channel flow simulations. When the 
computation for the pressure driven flow in a square duct was carried out using the D3Q19 formulation, a 
similar stability bound was encountered. 

Ideally, one would like to use a fixed value of r for the entire range of 0 < A < 1 in a simulation. 
Computational stability would then require the use of r around 0.6, instead of a value that is close to 0.5, 
which makes it difficult to simulate a lower viscosity, or higher Reynolds number flow. To overcome the 
restriction imposed by the numerical stability requirement due to interpolation, it would be useful if one 
could decrease the value of \ — (2 A — 1 )/r given by Eq. (1.11a). This can be accomplished by using 

(2-4) u bf = (^1-— )«; + — and X = for A > 1/2. 

That is, the velocity u^j is evaluated at ( x & + e Q /2), instead of at using the information at xj and x w 

through linear extrapolation. 

With Eq. (2.4) replacing Eq. (1.11a), the channel flow simulations using D3Q19 lattice model are carried 
out again for A from 0.85 to 1.0. Satisfactory results for the velocity profiles are obtained for r = 0.505 with 
N x x N y x N z = 5 x 35 x 5 in terms of computational stability. For A < 0.85, the accuracy of the solutions 
using Eqs. (l.lla) and (2.4) is the same when the computations are stable. 

3. Results and Discussions. 


3.1. Fully Developed Flow in a Square Duct. For fully developed flow inside a square duct of 
height H defined by the region —a<y<a, and— a < z < a. where a = H / 2, the axial velocity profile can 
be found in Ref. [28, p. 123]: 


(3-1) 


u x {y, z) = 



e»h (2^) 1 

cosh _ 


cosh (etgi£.) 

(2*+ 1)3 


Figure 3 compares the exact axial velocity profiles at 2 = 0 and the LBE-based solution with A = 0.2 and 
H = 2a = 32.4. A total of N x x N y x N~ = 13 x 35 x 35 grid points are used. The pressure gradient is 
^ — —1.0 x 10 -6 and r — 0.52. The nineteen- velocity model is used in the simulations. Excellent agreement 
was obtained. 

Figure 4a shows the dependence of relative i^-norm error, 


(3.2) 


E 2 = 


\ 


n n 

[uhBE{y, z) - U exact (y, z)] z dydz 

_ 

n H ’ 

[u ex!U ;t(y, Z)f dydz 

_ 


on the duct height or the lattice resolution H = N y — 3 4- 2 A. The integral is evaluated by the trapezoidal 
rule. As was demonstrated by Mei et al [20], the boundary treatment results in second order convergence 
for 2-D channel flow. Fig. 4a clearly shows that the total error (from both the flow field and the boundary 
condition) of the LBE solution in 3-D flow decays quadratically. 

Figure 4b shows the relative L 2 -norm error E 2 as a function of A in the duct flow using 13 x 35 x 35 
grid points and r = 0.52. For the purpose of comparison, the relative L 2 -norm error in the 2-D channel flow 
simulation using the D2Q9 model with N y — 35 and r = 0.52 is also shown. The relative error is larger in 
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u x(y) 

Fig. 3. Comparison of axial velocity profiles in a pres sure- driven square duck flow at z — 0 between the exact solution 
( solid line) and the LBE solution ( dashed line) with A = 0.2, dp/dx = — 1.0 ■ 10 -f \ r = 0.52, and H — 32.4. 

3-D duct flow than in the 2-D channel flow. Nevertheless, the error exhibits the same qualitative behavior 
in both 2-D and 3-D as a function of A. 

It should be noted that the accuracy of the D2Q9 model and the D3Q19 model is different in the sense 
that beyond the conserved moments (density and momentum in athermal fluids), these two models have 
different accuracy in preserving higher order moments (fluxes) [13, 14]. The D2Q9 model preserves all the 
moments up to the second order in momentum space, which include momentum fluxes, and maintains the 
isotropy of these moments, whereas the D3Q19 model can preserve density and momentum, but cannot 
maintain the same accuracy and isotropy of the fluxes as the D2Q9 model. The only 3-D equivalent of the 
D2Q9 model in terms of accuracy of the moments is the D3Q27 model [13, 14]. 

3.2. Simulation Results for 3-D Lid-Driven Cavity Flows. Lid-driven cavity flow has been stud- 
ied extensively in the CFD community. Most research has been focused on 2-D problems. Limited numbers 
of reliable numerical results for steady state 3-D cavity flows have been obtained in the past several years. 
In this study, the multi-block finite difference solution of the NS equations obtained recent by Salom [26] is 
used to compare with the present LBE based results. 

The size of the cavity is iJ 3 , the number of grid is N x x N y x N z , and N x = N y = N z . The driving lid 
is placed at y = H, moving along the direction of x-axis with a speed U = 0.1 in lattice units. Figure 5a 
compares profiles of horizontal velocity u x (y) obtained using 33 x 33 x 33 lattices with the solution to the 
NS equations at x/H — z/H — 0.5 for Re = 400. All three LBE models (D3Q15, D3Q19, and D3Q27) 
are used. For the fifteen-velocity model, the computation becomes unstable and blow's up at this Reynolds 
number with 33 3 lattice resolution and A = 0.5. For A = 0.5, the nineteen-velocity model and the twenty- 
seven- velocity model give very similar u x (y) profiles and both under-predict slightly the magnitude of the 
minimum in the profiles. The 19- velocity model is also used with A = 0.25; there is a slight overshoot in the 
velocity profiles in comparison to the results in Ref. [26]. Fig. 5b compares u x (y) profiles obtained using the 
fifteen-velocity and 19-velocity lattice models on the 67 3 lattice grids and A = 0.5 with the NS solution [26] 
at x/H = z/H = 0.5 for Re = 400. Excellent agreement is observed. Clearly, the nineteen-velocity model is 
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A 

Fig. 4. Error behaviors in a pres sure -driven square duck flow with dp/dx = —1.0 ■ 10 -6 , r = 0.52. (a) Dependence of the 
relative Li-norm error on the lattice resolution H = N y - 3 + 2 A, with A = 0.2. The straight line is a least-square fit of the 
data ( symbols ) with a slope of —2. (b) Relative L, 2 -norm error as a function of A for the 3-D duct flow ( solid line) and the 
2-D channel flow simulations . 


superior to the fifteen-velocity model. Although the fifteen-velocity model requires 21% less CPU time and 
storage than the nineteen-velocity model per lattice, it is not as robust as the nineteen-velocity model and 
may actually require more CPU time and memory to obtain a reasonable solution since more lattice points 
are clearly needed. 

It should be noted that stability property of the nineteen-velocity model and the fifteen-velocity model 
are significantly different. All LBE models have inherent spurious invariants because of their simple dynamics 
[18]. However, the stability of the LBE models, which is very much affected by these spurious invariants, 
differs from one model to another, and also depends on other factors such as boundary conditions and the 
local Reynolds number [18]. Among the three 3-D LBE models (D3Q15, D3Q19, and D3Q27), the D3Q15 
model is the least isotropic and therefore is more prone to numerical instability. This is independently 
verified in a recent work by Kandhai et ah [16]. It was observed that the D3Q15 model may induce artificial 
checkerboard invariants which are the eigenmodes of the linearized LBGK collision operator at wave vector 
k = 7r; this can cause spatial oscillations to develop in the flow field at high Reynolds number [18]. Although 
it was pointed out that the presence of solid walls can suppress the oscillation in certain cases, the solid 
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u x {y)/u 



u x {y)/u 

Fig. 5. 3-D lid-driven cavity flow. Normalized horizontal velocity profiles u x (x, y, z)/U at x/H = z/H = 0.5, with 
rm Re = 400, U = 0.1. The Namer-Stokes solution (solid line) from Ref. [26 j are compared with the following LBE solutions. 
(a) N x x Ny x N z = 33 3 , the nineteen-velocity LBE solutions with A = 0.5 ( dashed line ) and A = 0.25 (long dashed line), and 
the twenty-seven-velocity LBE solutions with A = 0.5 { dotted line), (b) N x x N y x N z = 67 3 , the nineteen-velocity ( dashed 
line ) and fifteen-velocity ( dotted line) LBE solutions with A = 0.5. 


walls in the present case actually excite the oscillation by producing a shear stress singularity at the two 
corners between the moving and stationary walls. Clearly, the D3Q19 model is better suited to handle flow 
singularities than the D3Q15 model in this case. 

Figure 6a compares the profiles of transversal velocity u y (x) obtained from various 3-D LBE models 
using 33 3 lattices (grids) with the NS solution at y/H = z/H = 0.5 for Re = 400. For A = 0.5 we found 
that the results from the twenty-seven-velocity model deviate more from the NS results of Ref. [26] than 
the results of the nineteen- velocity model with the same resolution in the spatial region 0.1 < x/H < 0.6. 
Both models under-predict the extrema of the velocity profile compared to the NS solution of Ref. [26]. For 
A = 0.25, the results of the nineteen-velocity model slightly over-predict the extrema, also shown in Fig. 6a. 
However the difference is relatively smaller in both cases. Figure 6b show's similar results of velocity profile 
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with a resolution of 67 3 grid points and the same Reynolds number Re — 400. With 67 3 lattice resolution, 
the result of the fifteen-velocity model significantly differs from the result of the nineteen-velocity model and 
that of the NS solution in Ref. [26]. These comparisons further suggest that the nineteen- velocity model is 
better than the fifteen-velocity model in terms of accuracy and stability, and better than the twenty-seven- 
velocity in terms of computational efficiency. The nineteen- velocity model represents a good compromise in 
terms of both computational efficiency and reliability. 



x/H 



Fig. 6. 3-D lid-driven cavity flow. Normalized transversal velocity profiles u y (x , y , z)/U at y/H ~ z/H — 0.5, with 
rm He = 400, U = 0.1. The Navier-Stokes solution (solid line) from Ref. [26] are compared with the following LBE solutions. 
(a) N x x Ay x N z = 33 3 , the nineteen- velocity LBE solutions with A = 0.25 ( long dashed line) and A = 0.5 (dotted line), and 
the twenty-seven-velocity LBE solution with A = 0.5 ( dashed line), (b) N x x Ay x N z = 67 3 , the nineteen-velocity ( dashed 
line) and fifteen-velocity ( dotted line) LBE solutions with A = 0.5. 

Figures 7 and 8 show the effect of Reynolds number (from 100 to 2000) on the profiles of horizontal 
velocity u x (y) and transversal velocity u y {x) at x/H = z/H = 0.5 based on the D3Q19 model. For Re — 100, 
400 and 1000, A = 0.5 is used. It is worth noting that for Re = 2000, the system size of 67 3 , U — 0.1, and 
r = 0.50325, the LBE simulation with A = 0.5 eventually becomes unstable, although the steady state result 
of Re = 1000 is used as the initial condition for Re = 2000. When A = 0.25 is used on the 67 3 lattice system, 
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no computational instability occurs and the steady state solution is obtained. Weak spatial oscillation in the 
u x {y) and u y {x) velocity profiles was observed for Re = 2000, which indicates that further increase in Re 
would require better spatial resolution. It is also worth pointing out that when FH’s boundary condition [8] 
is used for Re = 2000 with A = 0.25, the solution eventually blows up even when converged results (based 
on the present boundary condition for A = 0.25) at Re = 2000 are used as the initial condition. 



Fig. 7. 3-D lid-driven cavity flow. Effect of Reynolds number Re on the normalized horizontal velocity profiles 

u x (x.y,z)/U at x/H = z/H = 0.5. The results are obtain by using D3Q19 LBE model with V = 0.1, Re = 100, 
N x x Ny x N z = 33 s , A = 0.5 ( solid line); U = 0.1, Re = 400, N x x N y x N z = 67 3 , A = 0.5 ( dotted line); U = 0.1, 
Re = 1000, N x x Ny x N z = 67 3 , A = 0.5 ( dashed line); and U = 0.1, Re = 2000, N x x N y x N : = 67 3 , A = 0.25 ( long dashed 
line). 



x/H 

Fig. 8. 3-D lid-driven cavity flow. Effect of Reynolds number Re on the normalized transversal velocity profiles 

u y (x, y, z)/U at y/H = z/H = 0.5. The results are obtain by using D3Q19 LBE model with l = 0.1, Re = 100, 
Nx x N y x Nz = 33 3 , A = 0.5 (solid line); U = 0.1, Re = 400, N x x N y x N z = 67 3 , A = 0.5 ( dotted line); V = 0.1, 
Re = 1000, N x xN y x N : = 67 3 , A =: 0.5 ( dashed line); and V = 0.1, Re = 2000, N x x N y x N- = 67 3 , A = 0.5 (long dashed 
line). 
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3.3. Fully Developed Flows inside a Circular Pipe. Figure 9 shows the discretized domain and 
the boundary nodes Xb (denoted by solid symbols) for flow inside a circular pipe of radius R = 9.5 lattice 
units. Geometrically, the LBE simulation of the pipe flow differs from that of the duct flow in that the 
fraction of the intersected link A is not constant over the entire boundary. As seen in Fig. 4b, computational 
error can vary with A in the duct flow and the difference in the error can easily be as large as a factor of four 
for 0 < A < 1. Furthermore, the error is the smallest when A is between 0.3 to 0.6. Hence, it is reasonable 
to expect that, the overall error in the solution will depend on the distribution of A in the entire set of A. 



Fig. 9. Schematic for the boundary nodes x 5 (solid symbols) in yz plane 


in the pressure-driven in a pipe of radius R = 9.5. 


Figure 10 shows the relative L^-norm error for the axial velocity profile defined as 


(3.3) 


Eo = 


\ 


^ ^ [^LBE(^ji %k) 'tiexa.cti'yj i )] 
{yj>Zk)€V 

^ ^ [^exact (i/jn 

(y>,2fc)en 


where Q is the set of the discrete lattice grids inside the pipe, as a function of radius R for R = 3.5, 4.5, 
5.5, 9.5, 13.5, 18.5 and 23.5. The pressure gradient is = -1.0 x 10 -6 and r = 0.52. It is noted that 
each simple summation in Eq. (3.3) is slightly less than the exact integration over the entire circle due to 
the discretization. To ensure that such a treatment does not affect the qualitative behavior of the error 
measurement, the centerline axial velocity, u c , is also compared with the exact solution and the error is 
defined as: 


(3.4) 


E c = 


i^c,LBE exact! 


K, 


exact | 


It is seen that E c behaves very similarly to E 2 and both are non-monotonic. This oscillatory behavior could 
be due to the difference in the distribution of A, which in turn results in the difference of the dissipation 
due to the interpolation around the boundary. Shown also in Fig. 10 is the error E 2 of the square duct flow 
solution (with A = 0.2) as a function of equivalent radius H /tt */ 2 , which exhibits a quadratic convergence. 
Despite the non-monotonic behavior, it can still be seen that on average, E 2 and E c decay quadratically 
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with increasing radius and the accuracy in the pipe flow simulation is comparable to that in the square duct 
flow simulation. 



R or H/ 7T 


Fig. 10. Variation of relative errors E 2 and E c in velocity profiles as a function of the pipe radius R or the equivalent 
radius H/ tt 1 / 2 for a square duct. Shown are the E 2 (R) (O) and E C (R ) (O) for the pres sure- driven pipe flow, and E 2 (H/n l/2 ) 
(+) for the pressure-driven square-duct flow. The straight line is a least-square fit of E 2 (H/k 1/2 ) with a slope of -2. 

Figure 11 shows the axial velocity profiles in the pipe for R = 3.5, 5.5, 9.5, and 13.5 in comparison with 
the exact solution. Even for a very small radius R = 3.5, the LBE solution agrees with the exact solution 
remarkably well. A noticeable discrepancy in the velocity profile at R - 9.5 is also observed in E 2 and E c 
shown in Fig. 10. 



Fig. 11. Comparison of the axial velocity profiles between the LBE solutions ( dashed lines) and the analytic solutions 
(solid lines) for the pressure- driven pipe flow. Shown in the figure (from left to right) are solutions with R — 3.5, 5.5, 9.5 , and 
13.5. dp/dx = -1.0 • 10" 6 , t = 0.52 


3.4. Simulation Results for a Uniform Flow over a Sphere. The conventional LBE scheme uses 
uniform meshes. Without local mesh refinement, it is difficult to compute the external flow over a blunt 
body efficiently since a large number of grid points in the far field will be wasted. As a first attempt, the 
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flow over a sphere is computed within a finite region in the transversal directions. 

As shown in Fig. 12, the outer boundary is placed at y = ±H/2 and z — ±H/ 2. At y = —Hj 2, the 
lattice is j = 2. The boundary conditions at j = 1 for / rt \s are given by the following linear extrapolation: 

(3-5) /a(i\ 1, fc) = 2/ Q (i, 2, *) - / Q (i, 3, Ar). 

The velocity at j = 2 is set as 

(3.6) u(i, 2, fc) = u(i, 3, fc). 

Similar treatment is applied at = +H/2 and z — ±H/2. The extrapolation condition given by Eqs. (3.5) 
and (3.6) allow the flow to leave the outer boundary. This helps to reduce the effect of the outer boundary 
on the flow field and on the drag force. At the inlet, a uniform velocity profile is imposed at i — 1.5 (half way 
between the first and second lattice points) and Eq. (1.9) is applied to obtain the condition for / Q (l, j, k) 
with \ = 0. At the exit, a simple extrapolation is used, 

(3-7) f a (N x , j , k) = 2f a (N x - 1, j : k) - f a (N x - 2, j, k). 

On the surface of the sphere, Eqs. (1.9), (1.10), (1.11b), and (2.4) proposed in this work are used to update 
the boundary conditions for / Q ’s. Only the LBE nineteen- velocity model is used to simulate the flow over a 
sphere. 



Fig. 12. Schematic for computational domain in the uniform flow past a sphere. 


Figure 13 shows the velocity profile u x (y) based on a series of computations carried out for several values 
of the radius R = 3.0, 3.2, 3.4, 3.6, 3.8 and 4.0 for H/R = 10 at Re = 10. The results are obtained with 
r = 0.7. Fig. 14 compares the axial velocity profile (at y = z = 0) for the same set of parameters. It is 
worth noting that the present LBE computation does not have sufficient resolution for the given Reynolds 
number. Yet the velocity profiles agree with each other accurately. The fact that we have obtained a spatially 
accurate solution over a range of radii strongly suggests that the present boundary condition treatment for 
curved geometry in the LBE method is capable of handling more complex geometries while maintaining good 
accuracy. 

4. Concluding Remarks. Three 3-D LBE models, including the fifteen- velocity, the nineteen-velocity, 
and the twenty-seven-velocity model, have been assessed in terms of efficiency, accuracy, and robustness in 
lid-driven cavity flow. While accurate 3-D results can be obtained by using various LBE models, the nineteen- 
velocity model is found to be the best for the cases investigated. The fifteen- velocity model exhibits velocity 
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Fig. 13. Uniform flow over a sphere at Re =10 , r = 0.7, and H/R = 10. Comparison of the normalized velocity profiles 
u T {y)/U at x = z = 0 for various values for the sphere radius R = 3.0 (x), 3.2 (A), 3.4 (+), 3.6 (O), 3.8 (dashed line), and 
4.0 ( solid line). 



x/R 


Fig. 14. Uniform flow over a sphere at Re =10 , r = 0.7, and H/R = 10. Comparison of the normalized centerline at 
velocity profiles u x {x)/U (at y = z = 0) for various values for the sphere radius R = 3.0 ( dash-dot line), 3.2 (dotted line), 3.4 
( dashed line), 3.6 ( dash dot dot dot line), 3.8 (long dashed line), and 4.0 (solid line). 

oscillations and is prone to computational instability. The more complicated twenty-seven- velocity model 
does not necessarily give more accurate results than the nineteen-velocity model with the same spatial 
resolution. In this study, we have also modified the boundary condition treatment for LBE method proposed 
by Filippova & Hanel [8] and Mei et al. [20] when the fraction of the intersected link on the boundary A is 
greater than one half. This improves the computational stability when A is close to 1 and r close to 1/2. 

The simulations for flows in a square duct and in a circular pipe indicate that the current boundary 
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condition treatment for curved geometries results in second-order accuracy in 3-D flows. The velocity profiles 
for flow over a sphere show good self-consistency of the solution over a range of sphere radii used. 
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Appendix A. 6-U LBh Models. 

The D3Q15 model has the following set of discrete velocities: 

( (0> 0), a = 0, rest particle, 

(±1, 0, 0)c, (0, ±1, 0)c, (0, 0, ±l)c, a = 1, 2, . . . , 6, group I, 

[ (±1, ±1, ±l)c, a = 7, 8, . . . , 14, group III, 


Cq = < 


(AT) 

and the weighting factor w a is [24]: 





r 2/9, 

a = 0, 

rest particle, 


(A. 2) 


w a = 

] 1/9, 

«= 1, 2,..., 

6, group I, 





{ 1/72, 

a = 7, 8, . . . , 

14, group III. 


The D3Q19 model has the following set of discrete velocities: 




( (0,0,0), 



a = 0, 

rest particle, 

(A. 3) 

e a = < 

(±1, 0, 0)c, (0, ±1, 0)c, (0, 0, ±l)c, 

d — 1> 2, . . . , 6, 

group I, 



[ (±1, ±1, 0)c, (0, ±1, ±l)c, 

(±1, 0, ±l)c, 

a = 7, 8,..., 18, 

group II, 

and the weighting factor w a is [14]: 







r 1/3, 

a = 0, 

rest particle, 


(A. 4) 


w a = < 

1/18, 

a = 1, 2,.... 

> 6, group I, 





, 1/36, 

a = 7, 8, — 

, 18, group II. 


The D3Q27 model has the following discrete velocities: 





' (0,0,0), 



a = 0, 

rest particle, 

(A- 5) 

e a = < 

| (±1, 0, 0)c, (0, ±1, 0)c, (0, 0, ±l)c, 

a = 1, 2,..., 6, 

group I, 


(±1, ±1, 0)c, (0, ±1, ±l)c, 

(±1, 0, ±l)c, 

a — 7, 8, . . . , 18, 

group II, 



(±1, ±1, ±1 )c, 



q= 19, 20, , 26, 

group III, 

and the 

weighting factor w a is [14]: 






[ 

8/27, 

a = 0, 

rest particle, 


(A.6) 


w a - < 

2/27, 

Q = 1, 2,..., 

6, group I, 




1/54, 

a = 7, 8, . . . , 

18, group II, 




1 

1/216, 

a = 19, 20,.. 

. , 26, group III . 



In the above, c = SxjSt , Sx and St are the lattice constant and the 
structures for the D3Q15 and D3Q19 models are shown in Fig. 15. 
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Fig. 15. Discrete velocity sets {©a} for the D3Q15 (left) and D3Q19 (right) models. 
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